Se contruyo la base de datos con tres campos: fecha, tasa representativa del mercado y precio del barril del petróleo tipo brent, dado que este desde 2008 es el precio de referencia para la compañia Ecopetrol. La base de datos se aumenta utilizando variables indicadoras para los dÃas de la semana, los meses del año. La fuente de información para la serie del Brent es a través de U.S. Energy Information Administration (EIA), y para la tasa representativa del mercado (TRM) es el Banco de la República (BanRep) cargamos estas series a las cuales le asignaremos el nombre de Base.
#install.packages("quantmod")
#install.packages("tidyverse")
#install.packages("dplyr")
#install.packages("ggplot2")
library(quantmod)
library(tidyverse)
library(ggplot2)
library(dplyr)
library(readr)
library(lubridate)
library(stats4)
library(fUnitRoots)
library(aTSA)
library(timeSeries)
library(timeDate)
library(timsac)
getwd()
[1] "C:/Users/juanhena/Documents/SURA"
#setwd("/Users/jfhenaoduque/Google Drive/SURA")
setwd("C:/Users/juanhena/Documents/SURA")
# Base de datos
# PUNTO 1 -----------------------------------------------------------------
Base<-read.csv("BDPunto1.csv", sep=",", dec=".")
BDPunto1 <- read_delim("BDPunto1.csv", ",",
escape_double = FALSE, trim_ws = TRUE)
BDPunto1 <- separate(BDPunto1,"Fecha",c("day","month","year"), sep= "/")
#mensual
BDPunto1_m <- BDPunto1 %>% group_by(month,year) %>% summarize(mean_trm = mean(trm),
mean_brent = mean(brent, na.rm = T))
BDPunto1_m <- BDPunto1_m %>% arrange(year)
#anual
BDPunto1_y <- BDPunto1 %>% group_by(year) %>% summarize(mean_trm = mean(trm),
mean_brent = mean(brent,
na.rm = T))
En el código anterior, se entenderá que BDPunto1 es el nombre que tomará nuestra serie diaria. La base BDPunto1_m será nuestra serie en frecuencia mensual y la base BDPunto1_y será nuestra serie con frecuencia anual Como al momento de cargar los datos, Estos se verán como una hoja plana de información, es decir, el programa leerá los datos, pero asumirá que es una impresión, no manipulable; se contrarresta este inconveniente introducimendo el comando attach que nos ayudará a seccionar la información, dándole la opción al programa de transformar la información como si fuera una serie cargada en Excel, brindándonos completa manipulación.:
attach(Base)
#attach(BDPunto1)
#attach(BDPunto1_m)
#attach(BDPunto1_y)
Otro problema identificable en R es el hecho que por defecto los resultados del programa siempre seran vistos en formato de notación cientifica cuando los datos tienen muchos decimales. Entonces, para modificar los resultados, a fin de obtener el resultado con los decimales completos usamos la función:
options(scipen = 999)
Sobre la descarga de los datos: Inicialmente intente hacer la descarga del Brent por Yahoo Finance con frecuencia diaria del tickey BZ=F pero este no fue posible, ya que solo me generaba los últimos 20 datos, para testear si era algo del código o error mÃo (capa 8) lo probe con el del WTI que tiene etiqueta CL=F y la descarga si se efectuaba. No obstante, dado que la economÃa colombiana tiene como precio de referencia desde 2008 el Brent, realice su descarga en EIA que recoge las cifras arrojadas por Thomson Reuters
Realizamos los gráficos de la TRM y el precio del petróleo Brent en sus diferentes frecuencias:
plot(as.ts(BDPunto1$trm))
plot(as.ts(BDPunto1$brent))
plot(as.ts(BDPunto1_m$mean_trm))
plot(as.ts(BDPunto1_m$mean_brent))
plot(as.ts(BDPunto1_y$mean_trm))
plot(as.ts(BDPunto1_y$mean_brent))
Ahora gráficamos los datos con el fin de darnos una idea sobre los comportamientos de la TRM y el precio del brent library(plotly)
Base$Fecha<-as.Date(Base$Fecha,"%d/%m/%y")
library(plotly)
package 㤼㸱plotly㤼㸲 was built under R version 3.6.3Registered S3 method overwritten by 'data.table':
method from
print.data.table
Registered S3 methods overwritten by 'htmltools':
method from
print.html tools:rstudio
print.shiny.tag tools:rstudio
print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
Attaching package: 㤼㸱plotly㤼㸲
The following object is masked from 㤼㸱package:timeSeries㤼㸲:
filter
The following object is masked from 㤼㸱package:ggplot2㤼㸲:
last_plot
The following object is masked from 㤼㸱package:stats㤼㸲:
filter
The following object is masked from 㤼㸱package:graphics㤼㸲:
layout
plot_ly(data=subset(Base, subset= (Fecha>="2008-01-01")),
x = ~Fecha,
y = ~trm,
type= "scatter", mode="lines",
line=list(width=1,color="Blue"))%>%
layout(title="Tasa Representativa del Mercado",
xaxis=list(title="DÃa"),
yaxis=list(title="COP/USD"))
`arrange_()` is deprecated as of dplyr 0.7.0.
Please use `arrange()` instead.
See vignette('programming') for more help
[90mThis warning is displayed once every 8 hours.[39m
[90mCall `lifecycle::last_warnings()` to see where this warning was generated.[39m
Base$Fecha<-as.Date(Base$Fecha,"%d/%m/%y")
library(plotly)
plot_ly(data=subset(Base, subset= (Fecha>="2008-01-01")),
x = ~Fecha,
y = ~brent,
type= "scatter", mode="lines",
line=list(width=1,color="Red"))%>%
layout(title="Brent",
xaxis=list(title="DÃa"),
yaxis=list(title="Precio barril Brent (USD)"))
Utilizamos el diagrama de cajas y bigotes para explorar posibles relaciones
boxplot(BDPunto1_m$mean_trm)
boxplot(BDPunto1_m$mean_brent)
BDPunto1_m$month<-as.Date(Base$month,"%m")
plot_ly(data=subset(BDPunto1_m, subset= (BDPunto1_m$month<="01")),
x = ~month,
y = ~trm_mean,
type = "box")%>%
layout(title="TRM Mensual",
xaxis=list(title="Meses"),
yaxis=list(title="COP/USD"))
Con el fin de tener un análisis más completo incluimos los rezagos de la variable de interés
\[\begin{aligned} Corr(trm,brent)= Cov(trm,brent)/(S(trm)S(brent)) \end{aligned}\]
# Funciones -----
rezaga <- function(x,nlag){
rez <- cbind(as.ts(x),lag(x,nlag))
lagged <- rez[1:length(x),2]
return(lagged)
}
corsi <- function(y,x1){
correlacion <- matrix(data = NA, nrow = 1, ncol = 5)
correlacion[1,1] <- cor(y,rezaga(x1,0),use = "pairwise.complete.obs")
correlacion[1,2] <- cor(y,rezaga(x1,-1),use ="pairwise.complete.obs")
correlacion[1,3] <- cor(y,rezaga(x1,-2),use = "pairwise.complete.obs")
correlacion[1,4] <- cor(y,rezaga(x1,-3),use = "pairwise.complete.obs")
correlacion[1,5] <- cor(y,rezaga(x1,-4),use = "pairwise.complete.obs")
rownames(correlacion) <- c("Brent")
colnames(correlacion) <- c("0","-1","-2","-3","-4")
return(correlacion)
}
brent1<-as.ts(Base[-c(1:505),]$brent)
trm1<-as.ts(Base[-c(1:505),]$trm)
corsi(trm1,brent1)
0 -1 -2 -3 -4
Brent -0.8635655 -0.8637155 -0.8633608 -0.862887 -0.8623808
Dado que las dos series observadas sin importar su frecuencia presentan tendencia es pertinente realizarle pruebas de raÃces unitarias para determinar la pertinencia de hacer la correlación de estas dos variables niveles o si se deben realizar en retornos.
#install.packages("fUnitRoots")
adf.test(trm1)
Augmented Dickey-Fuller Test
alternative: stationary
Type 1: no drift no trend
lag ADF p.value
[1,] 0 1.71 0.978
[2,] 1 1.45 0.962
[3,] 2 1.41 0.959
[4,] 3 1.36 0.956
[5,] 4 1.43 0.961
[6,] 5 1.38 0.957
[7,] 6 1.42 0.960
[8,] 7 1.41 0.960
[9,] 8 1.41 0.960
Type 2: with drift no trend
lag ADF p.value
[1,] 0 0.2206 0.973
[2,] 1 -0.0172 0.954
[3,] 2 -0.0449 0.952
[4,] 3 -0.0982 0.946
[5,] 4 -0.0505 0.952
[6,] 5 -0.1089 0.945
[7,] 6 -0.0507 0.952
[8,] 7 -0.0557 0.951
[9,] 8 -0.0395 0.952
Type 3: with drift and trend
lag ADF p.value
[1,] 0 -2.40 0.407
[2,] 1 -2.64 0.305
[3,] 2 -2.70 0.283
[4,] 3 -2.76 0.258
[5,] 4 -2.66 0.298
[6,] 5 -2.71 0.277
[7,] 6 -2.68 0.291
[8,] 7 -2.69 0.285
[9,] 8 -2.70 0.281
----
Note: in fact, p.value = 0.01 means p.value <= 0.01
adf.test(brent1)
Augmented Dickey-Fuller Test
alternative: stationary
Type 1: no drift no trend
lag ADF p.value
[1,] 0 -0.776 0.401
[2,] 1 -0.763 0.406
[3,] 2 -0.730 0.418
[4,] 3 -0.720 0.421
[5,] 4 -0.736 0.416
[6,] 5 -0.736 0.416
[7,] 6 -0.767 0.405
[8,] 7 -0.748 0.412
[9,] 8 -0.763 0.406
Type 2: with drift no trend
lag ADF p.value
[1,] 0 -1.01 0.696
[2,] 1 -1.05 0.683
[3,] 2 -1.07 0.678
[4,] 3 -1.06 0.680
[5,] 4 -1.09 0.671
[6,] 5 -1.07 0.677
[7,] 6 -1.13 0.655
[8,] 7 -1.20 0.630
[9,] 8 -1.25 0.613
Type 3: with drift and trend
lag ADF p.value
[1,] 0 -2.20 0.494
[2,] 1 -2.26 0.466
[3,] 2 -2.34 0.434
[4,] 3 -2.35 0.428
[5,] 4 -2.35 0.429
[6,] 5 -2.33 0.437
[7,] 6 -2.34 0.432
[8,] 7 -2.46 0.383
[9,] 8 -2.49 0.370
----
Note: in fact, p.value = 0.01 means p.value <= 0.01
PP.test(trm1)
Phillips-Perron Unit Root Test
data: trm1
Dickey-Fuller = -2.6642, Truncation lag parameter = 9, p-value = 0.2971
PP.test(trm1)
Phillips-Perron Unit Root Test
data: trm1
Dickey-Fuller = -2.6642, Truncation lag parameter = 9, p-value = 0.2971
kpss.test(trm1)
KPSS Unit Root Test
alternative: nonstationary
Type 1: no drift no trend
lag stat p.value
12 0.0863 0.1
-----
Type 2: with drift no trend
lag stat p.value
12 0.119 0.1
-----
Type 1: with drift and trend
lag stat p.value
12 0.0473 0.1
-----------
Note: p.value = 0.01 means p.value <= 0.01
: p.value = 0.10 means p.value >= 0.10
kpss.test(brent1)
KPSS Unit Root Test
alternative: nonstationary
Type 1: no drift no trend
lag stat p.value
12 0.238 0.1
-----
Type 2: with drift no trend
lag stat p.value
12 0.273 0.1
-----
Type 1: with drift and trend
lag stat p.value
12 0.0929 0.1
-----------
Note: p.value = 0.01 means p.value <= 0.01
: p.value = 0.10 means p.value >= 0.10
Como se observa en las tablas, tanto para la prueba de Dickey-Fuller aumentada (ADF), Phillips-Perron (PP), y Kwiatowski, Phillips, Shy, Schmidt (KPSS) se encuentra que la serie tiene raÃces unitarias y por lo tanto se debe diferenciar la serie.
log <- data.frame(lapply(Base[,-1], log)) # saca logaritmos de determinado dataframe
names(log) <- paste("l",names(log),sep = "") # le da nombre a este nuevo dataframe
attach(log)
ret <- data.frame(lapply(log,diff,differences = 1)) # saca las diferencias de los logaritmos
names(ret) <- paste("r",names(ret),sep = "") # le da nombre a ese nuevo data frame
attach(ret)
Ahora tenemos la serie en retornos logaritmicos, con esto podemos procedemos a realizarle las pruebas de raÃces unitarias a la trm y al brent en rendimientos logaritmicos y ver si cumplen los supuestos de estacionariedad que nos permitan tener mayor precisión sobre la asociación lineal entre estas dos variables y potencialmente poder realizar estimaciones y hacer inferencia estadÃstica.
rlbrent1<-as.ts(ret[-c(1:504),]$rlbrent)
rltrm1<-as.ts(ret[-c(1:504),]$rltrm)
adf.test(rltrm1)
Augmented Dickey-Fuller Test
alternative: stationary
Type 1: no drift no trend
lag ADF p.value
[1,] 0 -45.3 0.01
[2,] 1 -35.1 0.01
[3,] 2 -28.9 0.01
[4,] 3 -25.9 0.01
[5,] 4 -22.7 0.01
[6,] 5 -21.3 0.01
[7,] 6 -19.8 0.01
[8,] 7 -18.3 0.01
[9,] 8 -16.5 0.01
Type 2: with drift no trend
lag ADF p.value
[1,] 0 -45.3 0.01
[2,] 1 -35.2 0.01
[3,] 2 -29.0 0.01
[4,] 3 -25.9 0.01
[5,] 4 -22.8 0.01
[6,] 5 -21.3 0.01
[7,] 6 -19.8 0.01
[8,] 7 -18.4 0.01
[9,] 8 -16.5 0.01
Type 3: with drift and trend
lag ADF p.value
[1,] 0 -45.4 0.01
[2,] 1 -35.2 0.01
[3,] 2 -29.0 0.01
[4,] 3 -26.0 0.01
[5,] 4 -22.8 0.01
[6,] 5 -21.3 0.01
[7,] 6 -19.9 0.01
[8,] 7 -18.4 0.01
[9,] 8 -16.6 0.01
----
Note: in fact, p.value = 0.01 means p.value <= 0.01
adf.test(rlbrent1)
Augmented Dickey-Fuller Test
alternative: stationary
Type 1: no drift no trend
lag ADF p.value
[1,] 0 -54.2 0.01
[2,] 1 -40.9 0.01
[3,] 2 -32.0 0.01
[4,] 3 -26.1 0.01
[5,] 4 -24.3 0.01
[6,] 5 -22.4 0.01
[7,] 6 -18.8 0.01
[8,] 7 -15.7 0.01
[9,] 8 -15.4 0.01
Type 2: with drift no trend
lag ADF p.value
[1,] 0 -54.2 0.01
[2,] 1 -40.9 0.01
[3,] 2 -32.0 0.01
[4,] 3 -26.1 0.01
[5,] 4 -24.3 0.01
[6,] 5 -22.4 0.01
[7,] 6 -18.8 0.01
[8,] 7 -15.7 0.01
[9,] 8 -15.4 0.01
Type 3: with drift and trend
lag ADF p.value
[1,] 0 -54.2 0.01
[2,] 1 -40.9 0.01
[3,] 2 -32.0 0.01
[4,] 3 -26.1 0.01
[5,] 4 -24.3 0.01
[6,] 5 -22.4 0.01
[7,] 6 -18.8 0.01
[8,] 7 -15.7 0.01
[9,] 8 -15.5 0.01
----
Note: in fact, p.value = 0.01 means p.value <= 0.01
PP.test(rltrm1)
Phillips-Perron Unit Root Test
data: rltrm1
Dickey-Fuller = -45.285, Truncation lag parameter = 9, p-value = 0.01
PP.test(rltrm1)
Phillips-Perron Unit Root Test
data: rltrm1
Dickey-Fuller = -45.285, Truncation lag parameter = 9, p-value = 0.01
kpss.test(rltrm1)
KPSS Unit Root Test
alternative: nonstationary
Type 1: no drift no trend
lag stat p.value
12 0.572 0.1
-----
Type 2: with drift no trend
lag stat p.value
12 0.163 0.1
-----
Type 1: with drift and trend
lag stat p.value
12 0.0751 0.1
-----------
Note: p.value = 0.01 means p.value <= 0.01
: p.value = 0.10 means p.value >= 0.10
kpss.test(rlbrent1)
KPSS Unit Root Test
alternative: nonstationary
Type 1: no drift no trend
lag stat p.value
12 0.0704 0.1
-----
Type 2: with drift no trend
lag stat p.value
12 0.071 0.1
-----
Type 1: with drift and trend
lag stat p.value
12 0.0466 0.1
-----------
Note: p.value = 0.01 means p.value <= 0.01
: p.value = 0.10 means p.value >= 0.10
Ahora con los retornos logarÃtmicos encontramos que el Brent y la TRM son estacionarias
corsi(rltrm1,rlbrent1)
0 -1 -2 -3 -4
Brent 0.01822319 -0.3350837 -0.08395169 -0.01031557 -0.03704103
Observese que ahora esta relación es positiva, y posiblemente bajo una estimación contemporanea el coeficiente estimado entre esas dos variables sea significativo. No obstante podemos pensar *** Ahora procedemos a obtener el promedio anual para la tasa representativa del mercado
attach(BDPunto1_y)
aggregate(mean_trm~year,data=BDPunto1_y, FUN=mean)
Ahora procedemos a obtener el promedio anual del precio internacional del petróleo
aggregate(mean_brent~year,data=BDPunto1_y, FUN=mean)
Ahora obtenemos el precio promedio para cada mes y cada año
aggregate(mean_trm~month*year,data=BDPunto1_m, FUN=mean)
aggregate(mean_brent~month*year,data=BDPunto1_m, FUN=mean)
BDPunto1_m$month<-strftime(BDPunto1_m$month, format = "%B")
BDPunto1_m$month<-strftime(BDPunto1_m$month, levels= month.name)
aggregate(mean_brent~*month,data=BDPunto1_m,FUN=mean)
Procedemos a representarlo en un gráfico
aggregate(mean_trm~year*month,data=BDPunto1_m, FUN=mean)%>%
plot_ly(
x = ~month,
y = ~mean_trm,
type = "scatter" ,mode = "lines",
split = ~year,
line=list(width=1))%>%
layout(title='Promedio diario mensual de la TRM',
xaxis=list(title="DÃas"),
yaxis=list(title="COP/USD"))
aggregate(mean_brent~year*month,data=BDPunto1_m, FUN=mean)%>%
plot_ly(
x = ~month,
y = ~mean_brent,
type = "scatter" ,mode = "lines",
split = ~year,
line=list(width=1))%>%
layout(title='Promedio diario mensual del Barril de Petróleo Brent',
xaxis=list(title="DÃas"),
yaxis=list(title="USD"))
La actividad económica colombiana ha sufrido notables cambios respecto a su materia prima de interés, casi durante todo el siglo XX el precio internacional de referencia para la economÃa era el precio del café, no obstante desde el descubrimiento de los yacimientos de Cusiana y Rubiales hizo que lentamente la actividad del petróleo tomará mayor relevancia para nuestra economÃa al constituirse en el principal producto de exportación. La importancia del precio del petróleo es tal que durante la crisis del petróleo en el 2014 se tuvieron consecuencias importantes para el paÃs , que obligaron a realizar reformas tributarias, ajustar el gasto estatal. Como se observa en la gráfica del Brent, este ha tenido caÃdas importantes en 2014 como consecuencia del descubrimiento en técnicas de Shale Gas y Shale Oil, y la sobreoferta de la OPEP. Durante este año, hubo notables caÃdas en el precio del petróleo producto de unas tensiones entre la OPEP en cabeza de Arabia Saudita y Rusia llevando a deterioros en la relación entre la OPEP y la OPEP+. No obstante, todo estas disputas se solucionarÃan en medio del confinamiento mundial producto de la pandemia del SARS-COV2, comunmente conocido como COVID-19.
Dada la fuerte importancia que tiene para Colombia las rentas petróleras y el flujo de dinero que llega de divisas como producto de la venta de barriles de petróleo el comportamiento de la tasa representativa del mercado (trm, en adelante) depende mucho del precio internacional. Si bien, hay una dependencia entre estas variableses bueno considerar si existe una posible cointegración de estos procesos, para evitar asà una posible relación espuria.
rel<-lm(ltrm~lbrent)
rel
Call:
lm(formula = ltrm ~ lbrent)
Coefficients:
(Intercept) lbrent
10.0478 -0.5316
Podemos señalar que de acuerdo a las pruebas de raÃces unitarias, y bajo lo propuesto por Granger en 1986 con las pruebas de cointegración se puede asumir que la trm y el precio del petróleo brent se encuentran cointegradas, es decir, hay una relación de largo plazo entre las dos, siendo está una elasticidad de largo plazo de -0.53. Aunque no es el objetivo del trabajo sugiero un modelo de mecanismo de corrección de errores. GRANGER, C.W.J (1986) ##"Developments in the Study of Co-Integrated Economic Variables"#. Oxford Bulletin of Economics and Statistics.
library(quantmod)
apple<-getSymbols("AAPL", from="2019-01-01", src="yahoo", auto.assign = F)[,6]
microsoft<-getSymbols("MSFT", from="2019-01-01", src="yahoo", auto.assign = F)[,6]
amazon<-getSymbols("AMZN", from="2019-01-01", src="yahoo", auto.assign = F)[,6]
facebook<-getSymbols("FB", from="2019-01-01", src="yahoo", auto.assign = F)[,6]
google<-getSymbols("GOOGL", from="2019-01-01", src="yahoo", auto.assign = F)[,6]
Datos.prueba<-data.frame(apple, microsoft, amazon, facebook, google)
Las series escogidas, fueron las tecnológicas, dada su alta valoración en el mercado y el interés que despiertan en especial por sus avances que han llevado a pensar que estas incursionen en los mercados financieros con sus propios bancos o hasta sus propias monedas como intento realizarlo Facebook con Libra, la cual ha tenido una serie de tropiezos desde su anuncio, en los que socios estrategicos del proyecto han decidido retirar su participación en la iniciativa de Facebook.
barChart(apple)
barChart(amazon)
barChart(facebook)
barChart(google)
barChart(microsoft)
En general observamos que todas estas series tuvieron una caÃda importante en el mes de marzo producto de la pandemia del COVID-19, aunque han tenido una notable recuperación, dado que este es uno de los sectores menos afectado por la pandemia e incluso algunas de estas empresas, Amazon por ejemplo, se ha visto favorecida por los confinamientos generalizados a nivel mundial debido a que su modelo de negocio se basa en las compras por internet. ## Apple
# Calculamos los rendimientos
applerends<-ROC(apple,type="discrete",na.pad = T)
applerendl<-ROC(apple,type="continuous",na.pad = T)
Datos.apple<-data.frame(apple,applerends,applerendl)
plot(Datos.apple$AAPL.Adjusted.1, type="l")
plot(Datos.apple$AAPL.Adjusted.2, type="l")
# Calculamos los rendimientos
amazonrends<-ROC(amazon,type="discrete",na.pad = T)
amazonrendl<-ROC(amazon,type="continuous",na.pad = T)
Datos.amazon<-data.frame(amazon,amazonrends,amazonrendl)
plot(Datos.amazon$AMZN.Adjusted.1, type="l")
plot(Datos.amazon$AMZN.Adjusted.2, type="l")
# Calculamos los rendimientos
facebookrends<-ROC(facebook,type="discrete",na.pad = T)
facebookrendl<-ROC(facebook,type="continuous",na.pad = T)
Datos.facebook<-data.frame(facebook,facebookrends,facebookrendl)
plot(Datos.facebook$FB.Adjusted.1, type="l")
plot(Datos.facebook$FB.Adjusted.2, type="l")
# Calculamos los rendimientos
googlerends<-ROC(google,type="discrete",na.pad = T)
googlerendl<-ROC(google,type="continuous",na.pad = T)
Datos.google<-data.frame(google,googlerends,googlerendl)
plot(Datos.google$GOOGL.Adjusted.1, type="l")
plot(Datos.google$GOOGL.Adjusted.2, type="l")
# Calculamos los rendimientos
microsoftrends<-ROC(microsoft,type="discrete",na.pad = T)
microsoftrendl<-ROC(microsoft,type="continuous",na.pad = T)
Datos.microsoft<-data.frame(microsoft,microsoftrends,microsoftrendl)
plot(Datos.microsoft$MSFT.Adjusted.1, type="l")
plot(Datos.microsoft$MSFT.Adjusted.2, type="l")
Se encuentra en otro archivo